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TURBULENT  BOUNDARY  LAYERS  OVER 
ROUGH  SURFACES  IN  HYPRSONIC  FLOW 

by 

John  M.  Russell 

Abstract 

A  method  for  predicting  the  downstream  development  of  momentum  thickness,  skin 
friction,  and  heat  transfer  in  a  supersonic  turbulent  boundary  layer  over  a  rough 
flat  plate  based  on  ideas  of  Van  Driest,  Rotta,  and  Bradshaw  is  derived  and  dis¬ 
cussed.  Admissable  thermal  boundary  conditions  include  the  case  of  prescribed 
wall  temperature  and  the  case  of  an  adiabatic  wall.  The  velocity  profiles  for 
compressible  nonadiabatic  flow  are  expressed  as  transformations  of  the  correspond¬ 
ing  velocity  profiles  in  incompressible  adiabatic  flow.  Analytical  curve  fits  to 
the  experimentally  determined  law-of-the-wall  (including  the  sublayer  region)  are 
given,  as  are  analytical  representations  of  the  effects  of  sand  grain  roughness 
based  on  the  well  known  data  of  Nikuradse. 

A  FORTRAN  source  code  for  implementing  the  method  is  included  as  are  sample 
calculations  of  the  momentum  thickness,  skin  friction,  and  heat  transfer  for 
several  roughness  heights. 
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I.  INTRODUCTION 

Interest  in  the  feasibility  and  technical  merit  of  a  high-altitude  high  speed 
cruise  missile  has  inspired  recent  Air  Force  interest  in  turbulent  boundary  layers 
over  rough  surfaces.  Apart  from  its  specific  military  applications,  the  problem 
is  of  fundamental  scientific  and  engineering  interest.  Though  the  subject  is  at 
least  as  old  as  the  supersonic  flight  of  aircraft,  the  physical  mechanisms  con¬ 
trolling  the  main  qualitative  features  of  such  flows  are  still  largely  a  matter  of 
debate,  and  there  is  ample  room  for  improvements  in  our  level  of  understanding  of 
the  physics  of  the  flows  and  the  equations  used  to  describe  them. 

A  long-range  experimental  program  under  the  direction  of  Dr.  Anthony  W.  Fiore  at 
the  High  Speed  Aero  Performance  Branch,  Wright  Aeronautical  Laboratories,  Dayton 
has  been  underway  for  over  two  years.  Among  the  goals  of  that  program  are  the 
aquisition  of  quality  experimental  data  contrasting  the  entropy  producing  mechan¬ 
isms  of  various  milled  surface  roughness  types  on  a  flat  plate  in  Mach  6  flow. 
Sometimes  the  interpretation  of  experimental  data  is  assisted  by  comparison  with 
numerical  predictions,  and  the  furnishing  of  such  predictions  was  the  original 
purpose  of  the  present  investigations. 

As  it  happens,  we  may  have  achieved  a  modest  simplification  and  improvment  of  the 
existing  technology  of  numerical  prediction  of  flat  plate  boundary  layers.  For 
example,  the  method  presented  herein  requires  only  that  an  explicit  first  order 
ordinary  differential  equation  be  solved  for  the  downstream  develoment  of  the  skin 
friction.  The  momentum  thickness  is  expressed  as  an  explicit  function  of  the  skin 
friction  (provided  the  parameters  of  the  free  stream  flow,  the  thermal  boundary 
conditions,  and  the  roughness  height  are  specified). 

II.  OBJECTIVES 

Our  objectives,  as  already  stated,  are 

(i)  To  provide  an  efficient  and  scientifically  plausible  algorithm  for  the  cal¬ 
culation  of  the  downstream  development  of  the  momentum  thickness,  skin  fric¬ 
tion  parameter  and  heat  transfer  parameter  for  realistic  input  data. 
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(ii)  To  develop  a  FORTRAN  source  code  to  implement  the  algorithm. 

(iii)  To  carry  out  sample  calculations  that  forecast  the  outcome  of  experiments 
underway  at  AFWAL 

III.  A  BOUNDARY  LAYER  CALCULATION  ALGORITHM 

In  this  section,  we  will  show,  starting  with  the  equations  of  the  motion  of  a 
visions  compressible  fluid,  that  the  downstream  development  of  a  turbulent  bound¬ 
ary  layer  over  a  rough  surface  in  supersonic  flow  can  be  modeled  effectively  by  a 
single  first  order  differential  equation  for  a  certain  friction  parameter.  The 
method  employs  a  few  well  known  results,  a  few  less  well  known  results,  and  a  few 
curve  fits  to  experimental  data  which  appear  here  for  the  first  time.  We  begin  by 
deriving  an  energy  equation  due  to  Rotta  (1959,  1960)  which  then  provides  the 
basis  for  a  compressible  law-of-the-wall  which  is  a  slight  generalization  of  the 
well  known  Van  Driest  transformation  (Van  Driest,  1951). 

A.  Rotta' s  Energy  Equation 

Let  (xj ^2X3)  denote  a  Cartesian  coordinate  system  with  x^  the  streamwise 
coordinate,  x2  the  vertical  coordinate  normal  to  the  wall  and  x^  the  spanwlse 
coordinate  defined  positive  in  the  right-handed  sense.  Let  (uj,  u2,  u^)  be  the 
local  instantaneous  velocity  field.  Let  p,  p,  and  T  be  the  instantaneous  pres¬ 
sure,  mass  density  and  temperature,  respectively.  Let  (q^,  q2,  q^)  be  the  local 
heat  flux  vector  due  to  thermal  conduction.  Let  p(T)  denote  the  molecular  vis¬ 
cosity  and  let  k(T)  denote  the  thermal  conductivity.  Let  0,  .  (j  =  1,2,3;  k  = 
1,2,3,)  denote  the  set  of  components  of  the  state  of  stress  tensor.  Let  Cp  and  cy 
denote  the  specific  heats  of  constant  pressure  and  at  constant  volume  respec¬ 
tively.  Let  R  *  Cp  -  cy.  Then  the  equation  of  state  for  an  ideal  gas  reads 

Pe  =*  PRT 

where  the  subscript  "e"  signifies  equilibrium  pressure  to  distinguish  it  from  the 
local  instantaneous  mechanical  pressure  p  defined  by 
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(Here,  as  elsewhere,  the  summation  convention  is  understood.)  These  two  types  of 
pressure  are  related  by 


P  -  Pe  =  - 


3u. 
PB  Ix^ 


where  u  is  the  bulk  viscosity.  Let 

D 


ekj 


1  A 

2  '  3x. 

J 


denote  the  elements  of  the  rate-of-def ormation  tensor.  Then  the  Stokes  viscosity 
formula  reads 


kj 


=  -n  6 


kj 


+  2u  [e 


kj 


mm 

3 


6.  .]. 

kj 


The  Fourier  heat  condition  law  reads 


Let  e  denote  the  thermometric  Internal  energy.  Let  R  denote  a  simply  connected 
region  in  space  and  6'  its  bounding  surface.  Let  dV  and  dE  denote  the  differential 
volume  element  and  the  differential  area  element  within  R  and  on  S ,  respectively. 
Finally,  let  (nj,  nj,  n^)  denote  the  components  of  the  outward  unit  normal  vector 
on  S.  Then  the  "control  volume"  forms  of  the  equations  of  conservation  of  mass 
momentum,  and  energy  are,  respectively. 


UJJpdV)  -  -  //  (P\)"k  d  £ 

R  S 


UH  p  uj  dvJ  * "  U  (puj)uknkdZ  +  ^0jkVz 


+  ^  Vkuj  dI  -  i/  Vk  dI 


provided  we  ignore  body  forces  such  as  gravity 
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6  )  denote  the  anisotropic  or  shear  stress  part  of  the  o,  . 

J  KJ  KJ 

matrix.  Then  in  view  of  previous  definitions,  the  Stokes  viscosity  formula  can  be 


Let  t  .  =  2p(e, . 


written 


3u 

a,  .  =  -  p  6  .  +  p  -r— 2-  S  +  t,  , 
ki  e  kj  B  3x  kj  kj 

m 


Let  the  ensemble  average  operator  be  denoted  by  pointed  brackets  <  >•  Assuming 

the  flow  is  statistically  stationary,  the  energy  equation  becomes 


0  =  ‘  ^<pUk(e  +  "f  +  "P")>nk  dE 


3u 

+  ii  <  <“b  IT  \j  *  Wnk  1,1  -  U<V"k dE  • 

o  m  S 

The  quantities  in  the  first  line  of  this  equation  represent  "bulk  mixing"  effects, 
namely  transport  of  internal  energy,  reversible  pressure  work  done  on  the  fluid  in 
the  control  volume  by  the  surroundings,  and  transport  of  kinetic  energy,  respec¬ 
tively. 


The  quantities  in  the  second  line,  by  contrast,  represent  "molecular  mixing"  ef¬ 
fects.  If  the  Reynolds  number  is  large  and  the  Prandtl  number  of  the  fluid  is  of 
order  unity  (as  it  is  for  air)  then  the  "molecular  mixing”  effects  will  be  impor¬ 
tant  only  in  the  sublayer  region  near  the  wall.  Under  these  same  hypotheses  the 
"bulk  mixing  effects"  will  be  dominant  in  the  interior  of  the  flow  away  from  the 
wall. 

AAA 

Let  (u,v,w)  =  (uj,u2,u3)  and  (x,y,z)  =  (xj^.x^).  Also  let  {i,j,k}  denote  the 
right  handed  orthonormal  triad  of  basis  vectors  associated  with  the  (x,y,z) 
coordinate  system.  We  now  consider  a  specific  choice  of  the  control  volume  H. 


Let  R  be  a  cross  sectional  fluid  slab  of  thickness  dx  in  the  x-direction,  of 
height  y  in  the  y  direction,  and  of  span  b  in  the  z-direction.  Assuming  air  is 
"calorically  perfect,"  we  have 
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Thus,  with 


e 


CPT* 


Then  on  the  top  face  of  the  slab  (where  n  =  j),  the  bulk  mixing  effects  dominate 
over  the  molecular  mixing  effects  so  that 


<pv(e  +  ~)>  =  c  <pvT>  »  -  <q0> 
p  P  %  2 


and 


u . u .  3u 

<P^V>  »  <UB  IT  62j  +  T2j)uj>- 

in 


We  will  assume  that  the  bottom  surface  of  the  slab  (n  =  -  j)  lies  on  the  wall  at 
the  level  y  =  0.  Also  the  wall  is  rigid  and  impermeable  to  mass.  We  also  assume 
that  the  flow  is  statistically  homogeneous  in  the  z-direction  so  that  the  net 

A  A  A  A 

energy  flux  across  the  faces  with  n  =  k  and  n  =  -  k  is  zero.  Then  the  energy 
equation  becomes 


u  u. 

0  =  -  Cp  <pvT>  -  <p  --L  jJ-v> 

+  ir  lly(-cp<euT>  -  <?u  +  <^B  IT  u> 

o  in 


+  <TMuj>  *  <91>)dy'] 


+  <q2> 


y-0* 


The  quantity  involving  the  x-derivative  represents  the  net  contribution  to  the 
energy  flux  across  the  forward  and  backward  faces  of  the  slab.  The  quantities  in 
the  first  and  last  line  represent  the  dominant  energy  flux  terms  across  the  top 
and  bottom  faces  of  the  slab,  respectively. 


In  the  "wall  region”  of  the  flow,  all  mean  flow  quantities  are  functions  of  the 
nondimens ional  variable  y+  defined  by 


Vw  ^visc 
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where  v  is  the  wall  value  of  the  kinematic  viscosity  and  u  is  the  friction 
w  T 

velocity.  For  boundary  layer  flows  in  general  and  for  flat  plate  flow  in  particu¬ 
lar,  i  visc  is  slowly  varying  in  x,  i.e. 


<n 


vise 


dx 


«  1 


We  will  accordingly,  neglect  the  x-derivative  terms  in  the  energy  equation  and 
obtain  the  leading  approximation. 

u.  u. 

-cp<pvT>  =  <p  -J^v>  +  qw 


where 


i  <q 


n>  n  =  -  <q0>  „ 

y=0  2  y=0 


The  next  step  in  deriving  Rotta's  energy  equation  is  to  apply  the  equations  of 
conservation  of  momentum  and  of  mass  in  order  to  rewrite  the  kinetic  energy  trans¬ 
port  term  in  a  more  useful  form. 


Assuming,  as  before,  that  the  flow  is  statistically  stationary  the  ensemble  aver¬ 
age  of  the  momentum  and  mass  equations  become 


du 


0  =  -  JJ  <pu.uk>nk  dl  -  //  <po>n.dE  +  //  <w 


B  3x 
m 


Tkj>nkdZ 


and 


respectively. 


0  =  -  //  <Puk>nk  dE  , 
S 


On  the  top  face  (n  *  j)  of 5  ,  the  "bulk  mixing"  terms  dominate  over  the  "molecular 
mixing"  terms,  so  that 


<pujV>  »  <V*B  if  {2j  +  V 
m 
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Assuming  the  flow  is  statistically  homogeneous  in  z  as  before,  the  momentum  flux 
balance  and  the  mass  flux  balance  for  a  cross  sectional  fluid  slab  become 


and 


0  =  -  <pujv>  -  <pe> <S2j 
3u 


+  *r  u0yt-<puju>  -  <pe>6  +  <wBl^  &  +  v]dy,] 

m 
3u 


+  <pe>y=0(S2j  “  <ViB  3x  62j  +  X2j>y=0 

m 


0  =  -  < pv>  -  —  [/rty  < pu>  dy’ J 


respectively.  We  assume,  as  before,  that  the  mean  flow  quantities  (except  the 
pressure)  are  functions  of  y+  =  y/ iv^sc  only  and  that 


vise 

dx 


«  1 


Also,  for  flat  plate  flow 


<pe>  =  P  =  constant. 


It  follows  that  the  x-derivative  terms  in  the  above  equations  are  zero.  Also, 
from  the  no-slip  boundary  condition  on  y  =  0  we  have 


<U 


— > 

B  3x  y-0 


3wv 

<yB  al>y=0 


0 


Letting 


<T21>y=0  "  V 


the  momentum  and  mass  conservation  equations  become 


0  -  <Pu.v>  +  <P£  |^>y=Q  62j  +  Tw  «Xj 


and 
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0  =  <pv>. 


respectively.  We  may  now  use  these  results  to  simplify  the  energy  equatior 


-c  <pvT>  =  <p  -*•  -l  v>  +  q 
P  2  ^w 


derived  earlier. 


Let  Uj  *  <uj> 

u  • '  =  u . -U • 
J  J  J 


u.u.  U.U. 


+  U.u.'  +  — rU- 
2  J  J  2 


so  that 


<P  v>  =  ~J2  < pv>  +  <pu^  ' v>  +  v> 

In  view  of  the  mass  conservation  equation  <pv>  =  0,  the  first  term  on  the  right 
vanishes.  Employing  both  mass  and  momentum,  the  second  term  becomes 

U.<pu.'v>  =  U-<pu.v>  -  U.U.<pv> 

J  J  J  3  j  J 

=  U.<pu.v> 

J  J 

=  -  V<n  |^>  -  UT 

B  3y  y=0  w 

where  U  =  Uj  and  V  =  Vj.  In  view  of  the  ordinary  boundary  layer  approximation,  V 
«  U,  so  the  second  term  on  the  right  dominates  over  the  first.  The  kinetic 
energy  transport  term  has  been  reduced  to 


U.u,  u.  u. 

<p  v>  =  -  Ui  +  <p  v> 

L  W  Z 


We  will  neglect  the  second  term  on  the  right  compared  to  the  first  since  we  expect 
third  order  products  of  fluctuations  to  be  small  compared  to  second  order  ones. 
With  this  approximation  the  full  energy  equation  becomes 
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-c_  <pvT>  =  -  Ut  +  q 
P  w  ^w 

which  we  will  call  Rotta's  energy  equation  (cf.  Rotta  (1959),  equation  31). 


B.  Van  Driest  Transformation  Theory 


We  now  introduce  the  "turbulent  Prandtl  number"  defined  by 


o 

t 


dT 

/ — <pu ' v ' >  \  r  dy  i 
1  dU~  n-<pvT>J 
dy 


where  T  =  <T>.  We  will  assume  that  is  a  constant.  Experimental  results 
suggest  =  0.9  is  a  sensible  average  according  to  Rotta  (1960). 


From  the  last  form  of  the  momentum  equation  written  down  in  the  last  section  (con¬ 
sidering  the  j  *  1  component)  we  have 


-<pu'v'>  *  T 


w 


(since  v'  *  v).  It  follows  from  the  above  definition  of  that 

-<>•'*>  ■  §  £ 

Rotta's  energy  equation  then  becomes 

cp  Tw  dT  „  .  . 

0C  dU  w  Mw 

Letting  pw  -  <P>y=Q 

TW  =  <T>y=Q 

T  »  P  U  2 
W  W  T 


d(f) 

w 

-(5-) 


2 

u.  a 


f  W  \  _  T 

u  T  c  '  °t  T  c 


t  U 


W  T  W  p 


U 

W  p  T 


we  have 
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Let 


0  i 


q  p  u  T  c 
w  t  w  p 


a  2  =  (c  -c  )T  -  c  (y-1)T 
w  C  p  V  w  p  w 

V  r 

u 


Then 


so  that 


Mi  — 
t  a 

2  2  2  V 
u  a  Mao  „ 

T  t  T  W  t  ..  2,  .  . 

T  c  f  2  “  Mx  (Y-1)ot 

w  p  [a  /  ( Y— 1 ) J 

w 


dg-) 


Vt  -  MT2(Y-DatC5-)  , 


hence 


M 


f  ■ c. +  V.  O  -  -f-  <’-»v^2 

f-T  •  T  J 


-  ^T1  MT2°t  {LC  - i?  +  ]  -  [J - ^ 

(Y-l)M/  H_LLm/o^  ut  ( Y~1 )  m  2 


CJ 

2  \  ~t 


rH 


It  Is  convenient  to  introduce  a  new  variable  $  defined  by 

0  0  „  C, 


U  _ q  _  [  (• _ q  1 2 _ 1  i  V?  .  . 

u  2  “  l ^  2'  Cy-1)  2  J  S^n  ^ 

t  (y-1)m/  (Y-1)m/  o. 


2  T  t 


Then,  in  terms  of  <j>. 


0  2o. 


(c.  +  — ^  ■■  j)  cos  $ 
2(y-1)M/ 


U) 


(2) 


The  above  formula  is  very  convenient  for  deriving  a  compressible  flow  generaliza¬ 
tion  of  the  logarithmic  form  of  the  law  of  the  wall. 

Let  6  denote  the  overall  thickness  of  a  turbulent  boundary  layer  (from  the  wall  to 
the  average  elevation  of  the  turbulent-non-turbulent  interface).  We  will  refer  to 
the  "generalized  logarithmic  region"  of  the  turbulent  boundary  layer  as  the  range 
of  values  of  y  which  satisfy,  say. 


(35)  lvLsc  <  y  <  (0.2)  6. 


Let  the  symbol  <  denote  the  Karman  constant.  Then  the  formula 


dU 

dy 


1_  V2 

*y  P(y) 


is  well  established  for  uniform  denlsty  flows.  Following  Van  Driest  (1951),  we 
will  assume  that  this  formula  also  holds  for  nonuniforro  mean  density. 


We  have  already  employed  the  hypothesis  that  the  mean  pressure  is  constant  (and 
uniform  in  space)  in  flat  plate  flow.  From  the  equation  of  state  for  an  ideal 
gas,  we  have,  therefore. 


p  *  (cp-cv)pT 

p  *  (cn"°»>P  T 
P  V  W  w 

(where  the  subscript  "w"  denotes  the  wall  level  (y»0).  It  follows  that 


T 

p  w 


an  identity  which  we  will  employ  frequently  in  what  follows.  With  the  definition 
-  2 

T  **  P  u  ,  the  formula  for  dU/dy  becomes 
w  w  t  J 


dU  _  UT(T(y)  s  l/2 
dy  ^y^  Tw  > 


or 


Our  original  definition  (1)  of  the  variable  $  above  was  a  formula  relating  U/u^ 
to  4>.  Rotta's  energy  equation  and  constancy  of  a  then  led  to  formula  (2)  which 
relates  T/T  to  It  follows  that  (1)  and  (2)  can  be  employed  to  eliminate  U/u 

w  T 
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and 

be 


T/T 

w 


from  the  right  hand  side  of  the  above  equation. 


I  iz  = 

<  y 


1 


d<t> 


The  result  is  found  to 


from  which 

4>  “  *0  -  Mt2  ot)  V2  (i  In  y+  +  C,,)  (3) 

(where  we  have  introduced  two  constants  of  integration  in  order  to  permit  one  to 
be  chosen  for  analytical  convenience  later). 


From  the  identity 


sin  <{>  =  sin  ( <p  —•  d>  +  0  ) 
o  o 


=  sin  C  4>  —  4>  )  cos  $  +  cos  ( <|>  —  <t>  )  sin  <p 

0  0  o  o 


Substituting  this  into  (1)  and  using  (3)  to  eliminate  4>  -  we  have 


x  (Y-1)M„ 


[( 


6. 


r)- 


'i 


V2 


•  (sin  [  Mt2  ot  )  1/2  (~  lny+  +  C2  )  ]  cos  <pQ 

+  cos  [  (— 2  Mt2  ot  J  1/2  (7  lny+  +  C2)J  sin  <f>o} 


(A) 


Since  both  of  the  constants  and  C2  arose  from  a  single  integration,  we  may 

choose  a  particular  value  for  either  one  without  loss  of  generality.  We  will 

choose  4>  such  that 
o 


l‘“.o 

q 

M  -►  0 


T 


+  c 


2 


which  will  ensure  that  our  formulas  (which  hold  for  compressible  flow  with  heat 
transfer)  reduce  to  the  appropriate  form  in  the  incompressible  adiabatic  limit. 
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If  we  write  (4)  in  the  form 


£_ 

u 


e  2  a 


[-  ^  C-2  *  C, )  V2 

2(Y-1)M/ 


•  MT2otr  1/2  sin  [  (— MT2°t)  1/2  (7  !ny+  +  C2)  ]  cos  ^ 


+  (iltil  MT2otrl/2cos  MT2ot)  1/2  (i  lny+  +  C2)]  sin  *q 


Pn  3  2af 

* — Hr  l— l-H*  C1)"1/2I 


(y-1)Mt"  2(y-1)Mt 

The  above  limit  condition  therefore  leads  to  the  choice 

&  V 


2(y-1)M 


sl"  *»  ■  MT2°t> 1/2 *  cil'‘/2 


e  V 


(y-Dmt"  2(y-i)mt 


2  2 

(which  incidentally,  satisfies  sin  4>q  +  cos  <f>Q  =  1) 


Eliminating  <j>o  from  the  above  formula  for  U/u^,  we  get 
3. 


l - 1<l-T  +  Cl1/2(-tol  MT2otr1/2sin  [{^1  MT2ot)1/2(i  lny+  +  C  )] 

ux  (Y-l)M  1  1  2  T  t  2  2 

3 

- 9 — _  cos 

(Y-1)m/ 


ll^  M,\) 1/2  (7  H  ♦  c2)] 


(cf  Rotta  1968,  equation  17). 


The  above  formula  relates  the  compressible  heat  conducting  velocity  profile  (in 
the  generalized  logarithmic  region)  to  the  incompressible  adiabatic  velocity  pro¬ 
file  (in  the  ordinary  logarithmic  region).  If  we  let 
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Utnc  <y+)  5  7  lny+  +  C2 

(in  the  ordinary  logarithmic  region),  then  the  above  compressible-incompressible 
transformation  formula  becomes 


T  (Y-l)M 


a-r*  1/2»L  <y+> 


(Y-l)M 


?TT~  l(^,2"t)1/2C  <y+» 


(5) 


C.  Von  Karman  Momentum  Integral 

The  third  general  result  which  we  will  employ  in  our  boundary  layer  calculation 
algorithm  is  the  Von  Karra/Jn  momentum  integral.  The  specific  form  which  we  wish  to 
employ  has  been  derived  by  many  authors  and,  unlike  the  Rotta  energy  equation  and 
compressibile-incompressible  transformation  formula,  we  feel  no  need  to  supply  a 
derivation  of  it  here.  Readers  are  referred,  instead,  to  the  derivation  of  Young 
(1953). 


Let  the  subscript  6  denote  conditions  at  the  outer  edge  of  the  boundary  layer.  We 
define  the  displacement  thickness  by 

s. 

pfi  s 

We  define  the  momentum  thickness  9  by 

0  =  J«  Sill  «<£>.(  j  -  Ujzl)  dy 
J°  U6  U6 

Then  the  Von  ICSrra^n  momentum  integral  is 


k  'W8’ 


T 

W 


with  a  variety  of  error  terms  among  which  are  the  apparent  stresses  at  the  outer 
edge  of  the  boundary  layer  due  to  accoustic  radiation  effects  and  the  effects  of 
x-derivatives  of  the  Reynolds  normal  stresses. 


D.  Consolodated  Algorithm  for  Smooth  and  Rough  Halls 

For  flat  plate  flow,  is  a  constant.  Inspection  of  equation  (6)  shows  that  if  0 
were  a  function  of  only,  then  one  could  eliminate  either  6  or  form  (6)  and 
obtain  a  first  order  ordinary  differential  equation  for  whichever  variable  was  not 
eliminated.  Although  it  is  hardly  apparent  at  this  stage,  we  will  show  in  the 
present  subsection  that  6  is  indeed  a  function  of  only  provided  representative 
data  is  supplied  that  would  normally  constitute  the  input  to  a  boundary  layer  cal¬ 
culation. 

(a)  Whole-Layer  Formula  for  ut  (y+,k+) 

inc 

We  begin  by  writing  down  two  analytical  curve  fits  to  well  known  experimental  data 
in  incompressible  flows  which  are  usually  presented  in  either  graphical  or  tabular 
form. 

Coles  (1956,  1968)  has  shown  that  nearly  all  turbulent  boundary  layers  (over 

smooth  walls)  can  be  fit  to  a  formula  of  the  form 

Uinc(y+)  mf(y+)  +  7  w({> 


(smooth) 
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in  which  f  and  w  are  ostensibly  universal  functions  determined  experimentally  and 
II  is  a  parameter  (called  "Coles'  wake  function").  A  tabulation  of  /(y+)  was  given 
by  Coles  (1955).  A  tabulation  of  w(y/6)  was  given  by  Coles  (1956).  An  analytical 
curve  fit  to  w(y/6)  was  given  by  Coles  (1968)  as 

w(f)  =  2  sin  2(-^)  =  1-cos  (-JP) 

In  the  "logarithmic  region,"  f( y+)  has  the  well  known  form 

/( y+)  =  7  lny+  +  c2  (7) 

In  the  sublayer  and  buffer  region  (say  y+  <  35)  the  above  formula  for  f( y+)  fails. 
The  true  law  of  the  wall  function  must  satisfy  the  constraint 


+ 

dy  y  =0 


0 


(8) 


Now  the  functions  tanh  (  )  and  sinh~^(  )  are  both  "concave  downward  for  all 
positive  argument  and  have  zero  concavity  at  the  origin — both  properties  of  the 
law  of  the  wall  function  /(y+).  Also,  any  linear  combination  of  sinh-1(  )  and 
tanh  (  )  will  be  asymptotic  to  a  logarithm  for  large  argument.  An  obvious  curve 
fit  to  ^(y+)  is  therefore  a  function  of  the  form 


tanh^) 


in  which  the  parameters  a,  d,  and  c  are  constrained  so  that  y+)  is  asymptotic  to 
the  form  (7)  for  large  y+  and  satisfies  the  slope  constraint  (8).  We  find  (for 
smooth  walls)  that 


d 


+ 


lna 

K 


d(1  -  27I>_1 


and  only  "a"  remains  adjustable.  By  trial  and  error,  we  have  found  that 


a  -  2.45 


-21- 


leads  to  a  uniformly  valid  approximation  to  Coles  (1955)  tabulation  (assuming  the 
values  of  <  and  C2  used  there)  with  a  maximum  error  under  two  percent. 

Combining  this  fit  to  (y+)  with  Coles  (1968)  fit  to  w(y/6),  a  whole-layer  formula 

for  ut  (y+)  (for  smooth  walls)  is  found  to  be 
inc 

ut  (y+)  =  —  sinh  *  ]  +  dtanh  (2—  J  +  —  ( 1-cos  (— ^-}J 

(smooth) 


with  the  above  values  of  a,b,  and  c. 


As  a  means  of  incorporating  roughness  effects,  we  will  apply  the  well-known  sand 
grain  roughness  data  of  Nikuradse  (cf  Cebeci  and  Bradshaw  (1977),  section  6.5). 
Let  k  now  denote  the  roughness  height  (not  to  be  confused  with  the  thermal  conduc¬ 
tivity  symbol  used  in  subsection  A  above).  Let 


k 


+ 


For  a  given  roughness  type,  the  velocity  profile  of  a  turbulent  boundary  layer 
exhibits  a  logarithmic  region  as  on  smooth  walls,  however  the  value  of  the  addi¬ 
tive  constant  C2  in  (7)  is  reduced  by  an  amount  Au+(k+).  An  analytical  curve  fit 
which  agrees  with  Nikuradse's  experimental  points  for  Au+(k+)  to  within  the 
scatter  of  those  data  is 


Au+(k+) 

with 


+  7  lnOi 


(k+)2 


(LsV  + 


+  2J 

(kV 


(9) 


(C2,  B2oo,  k,  Ls+)  =  (5.5,  8.5,  0. AO,  13.0) 

For  readers  wishing  to  test  this  formula,  we  note  that  the  Nikuradse  data  actually 
plotted  in  figure  6.16  of  Cebeci  and  Bradshaw  (1977)  is  the  function  B2(k+) 
defined  by 


B2(k+)  =  i  ln(k+)  +  C2  -  Au+(k+) 
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Si  nee  the  main  effect  of  sand  grain  roughness  in  the  logarithmic  part  of  the  velo¬ 
city  profile  is  to  reduce  the  additive  constant  C2  by  a  roughness  dependent  amount 
Au  (k  ),  it  seems  most  natural  to  incorporate  roughness  effects  into  our  whole- 
layer  velocity  profile  formula  by  an  identical  reduction  of  C2  at  the  point  where 
C2  appears  in  that  formula.  Specifically,  we  may  propose  a  whole-layer  velocity 
profile  formula  of  the  form 

U^nc(y+,k+)  =  ■■  sinh  1  +  d(k+)  tanh  (l-cos  (10a) 

where 

d(k+)  =  C2  -  Au+(k+)  +  —2.  (10b) 

and  Au+(k+)  is  given  by  (9)  above. 


(b)  Thermal  Boundary  Conditions 


In  the  subsection  titled  "Van  Driest  Transformation  Theory"  above  we  derived  the 
formula 
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f  -  ci*  v,(5-)  -  t-  ‘T-o  »t(fr 


Employing  the  subscript  S  to  denote  conditions  at  the  "outer  edge"  y  =  6  as  before 
and  employing  the  identity 


Hi  - M&a&  -  Hi  A  \ 

u  M  a  M  ^  T  ' 


u  M  a  M  '■T 

T  T  W  T  W 


we  have 


Hi  -  r  +  a  Hi  f!ii  V2  (Y-l)  _  M  2  Hi 

T  C1  +  Vt  M  I?  J  2  °tM6  T 


We  will  assume  that  the  given  data  constituting  the  input  to  the  boundary  layer 

calculation  will  include  a  thermal  boundary  condition  which  will  be  either  an 

adiabatic  wall  condition  3  =  0  or  a  given  wall  temperature  condition  T  =  given. 

q  w 

In  the  adiabatic  wall  case  (3^=0),  the  above  formula  yields  a  wall  temperature 
formula  of  the  form 


(TW)V°  =  C 


:(14¥  “A2* 


(11a) 


In  the  given  wall  temperature  case,  the  above  general  formula  may  be  solved  for 


0q  -  K  ft)"1  (rr1/2[^rU  +  °tMs2)  -  cj  (nw 

T  W  W 

All  of  the  parameters  in  (lib)  other  than  would  normally  be  known  before  the 
boundary  layer  calculation  is  started.  In  this  sense  (lib)  defines  a  function 
3  (M  )  for  given  T  .  If  £5  is  given  to  be  zero  then  (11a)  defines  T  .  In  any 

(|  T  w  (J  " 

case  the  quantities  3^  and  Tw  are  each  either  known  from  given  data  or  expressible 
as  an  explicit  function  of  M  . 


(c)  Determination  of  the  Functions  6(11^)  and  M^(x) 


In  formula  (1)  of  subsection  B  above,  we  introduced  a  variable  <f>  which  is  a 
function  of  U/u^.  Letting  denote  the  value  of  $  at  y  =  6  and  employing  the  now 
familiar  identity 
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u 

T 


M,  1/ 

M  ''T  ‘ 


w 


we  find  from  (1)  that 
-1 


sin 


U(- 


(Y-l)M. 


r)2  + 


2C . 


(Y-1)MT  oc 


|“  1/f2  [”6_  (Zi)  l/‘ 

Lm  '•T  1 


2  _ 


(Y-l)M, 


"]}  02) 


In  that  same  subsection  a  constant  of  integration  (which  first  appeared  in 

formula  (31))  arose  and  was  assigned  a  specific  value  as  a  result  of  the  require¬ 
ment  that  the  compressible  heat  conducting  form  of  the  law  of  the  wall  in  the 
generalized  logarithmic  region  reduce  to  the  known  incompressible  adiabatic  law  of 

the  wall  as  6  and  M  both  tend  to  zero.  Taking  the  arcsine  of  the  formula  for 
q  t 

sin<j>  so  determined,  we  have 
o 

<t>  =  sin  Ml( - 9 — j)2  + - ~i — 1  /2  1- - ■CL— 03) 

°  (y-Dmtz  (y-Dmt  ot  (y-Dm^ 

The  parameters  <j>g  and  4>o  are  now  known  functions  of  once  the  thermal  boundary 

conditions  and  the  normal  free  stream  data  of  the  boundary  layer  claculation  prob¬ 
lem  are  specified. 

Now  we  may  obtain  a  whole-layer  generalization  of  formula  (3)  of  subsection  B 
above  by  replacing  the  quanitity 


7  lny+  +  C2 

by  the  more  general  expression  U+nc(y+,k+)  as  defined  by  forumla  (10)  of  subsec¬ 
tion  D.  In  particular,  at  the  outer  edge  y  *  6,  we  get 

♦j  -  ♦„  *  1-2-11  "t2°t )  1/2  UL 

where,  from  (10), 

Ulnc  (6+,k+)  “  7  ln  6+  +  C2  -  *u+(k+)  +  ” 

Eliminating  uT  between  the  last  two  equations  and  solving  for  6+,  we  get 
Inc 
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6+  =  expl4(^-  MT2otr  V2  U,  -  *Q)  “  C2  +  Au+(k+)  --^j} 

In  view  of  previously  derived  results,  the  right  hand  side  is  a  known  function  of 
provided  the  usual  free  stream  data  and  thermal  boundary  conditions  for  a  flat 
plate  boundary  layer  calculation  problem  are  specified. 

The  function  0(M^)  is  now  almost  completely  determined.  The  only  remaining  func¬ 
tion  we  require  is  the  function  J  visc^t^  defined  by 


7  (M  )  =  — 

‘-Vise'  T  U 


*<V 

p  a  M 

W  W  T 


RT 


"<V 


p  RT  (yRT  ) 
w  w  w 


T/2m 


.  RT  ,,  p(T  ) 

l_  i _ Wn  l/2  w 

P  ^  Y  J  M 

i 


where  P  is  the  pressue  and  u(Tw)  is  the  wall  value  of  the  molecular  viscosity 
u(T).  From  the  "Southerland  law",  we  have 


b(T) 


ref  '•T 


3/2 


ref 


Tref  +  SI 
T  +  S, 


where  (for  air) 


W  =  (0.350)  10"6 
ref  ft2 

Tref  "  492  °R 
Sj  =  198  °R 

From  the  definition  of  the  momentum  thickness 

e  ■  j0s  m  (1  -  5^) 


and  the  identities 
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ri'i*  • 

W  T  T 

u  M  T  i.  . 

1  U  _  T  ,  w>  V?  U 

Ur  u  Mr  '■Tr'  u  * 


u_ 

U, 


and  y 


y  11  vise'  we  have 


6+(Mt) 


e(M  )  -  i  . 

T  V1SC 


<V  l  l£  (r) 1/2  lc! +  Vt  h  -  ~r  <’(-1)0t 


•  —  [1  -  ^  Or)''2  — ]1  dy+ 
u  L  M,  'T,'  u  J  J 


Given  the  type  of  input  data  necessary  for  boundary  layer  calculation  problems 
including  free  stream  conditions  and  thermal  boundary  conditions  and  given  the 
functions  of  already  defined,  the  right  hand  side  of  the  above  formula  for 
0(Mt)  is  fully  determined. 


In  particular  the  function  dB/dM^  is  a  known  function  of  Mt  (though  its  evaluation 
would  normally  require  numerical  differentiation).  It  follows  that  the  Von  KarmSn 
momentum  integral  (cf.  section  C  above) 


d_9  A-,2 

dx  -  ^ 

may  be  divided  by  dO/dM^  to  give 


which  is  a  first  order  autonomous  ordinary  differential  equation  for  as  a 
function  of  x.  It  is  only  necessary  to  specify  at  an  initial  value  of  x  to 
determine  uniquely  the  function  M^(x)  for  all  downstream  values  of  x. 


IV.  CALCULATION  EXAMPLE 


In  this  section,  we  will  apply  the  algorithm  developed  in  the  preceeding  section 
to  a  phsyical  problem  intended  to  model  a  boundary  layer  experiment  currently 
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underway  at  the  High  Speed  Aero  Performance  Branch  of  the  Air  Force  Wright 
Aeronautical  Laboratory  in  Dayton. 

A.  Physical  Setup 

Consider  a  supersonic  wind  tunnel  with  air  as  the  working  fluid.  Hie  tunnel  res¬ 
ervoir  conditions  are  denoted  by  the  subscript  ”o”.  TQ  and  PQ  are  both  given. 
From  the  design  of  the  tunnel  throat  and  test  section,  the  test  section  Mach  num¬ 
ber  is  controllable  and  is  therefore  taken  as  a  given  quantity.  Assuming  the 
flow  between  the  reservoir  and  the  test  section  is  isentropic,  we  have 


which  determines  the  "free  stream"  temperature  and  pressue  in  the  test  section  in 
terms  of  the  given  data. 

We  will  assume  that  the  model  consists  of  a  roughened  flat  plate  at  zero  inci¬ 
dence.  The  leading  edge  of  the  plate  coincides  with  the  straight  line  x=o,  y=o  in 
the  (x,y,z)  coordinate  system.  The  free  stream  velocity  is  in  the  direction  of 
the  positive  x  axis. 

The  plate  is  smooth  in  the  strip  0  <  x  <  immediately  behind  the  leading  edge. 
Behind  that  strip  the  plate  has  a  uniformly  distributed  roughness  of  height  k. 

B.  The  Starting  Laminar  Boundary  Layer 

The  development  of  a  flat  plate  laminar  boundary  layer  in  a  compressible  flow  can 
be  calculated  analytically  (cf.  Schlichting  (1968),  Chapter  XIII)  provided  the 
viscosity-temperature  relation  has  the  idealized  form 

»<T>  -  (rr)* 

ref 

where  u  and  Trej  are  constants.  Notice  that  the  above  formula  implies 
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U(Tfi) 

p(V 


(T./T  )“ 

o  w 


If  Tg  is  given  and  Tw  is  determined  from  the  thermal  boundary  conditions  of  the 
problem  then  w(T^)  and  y(Tw)  can  be  calculated  from  the  Southerland  law.  The 
above  formula  then  determines  w. 


For  example,  if  T0  =  1100  °R  and  =  6  then  from  the  isentropic  relations  = 

134.15  °R.  If  the  wall  is  adiabatic  and  the  "turbulent  Prandtl  number"  is  taken 

to  be  o  =  0.9  (Rotta  (1960),  then  from  the  formula  for  T.  in  section  D(b)  above, 
t 

we  get  Tw  =  1003  °R.  From  the  Southerland  law,  the  values 

p(T5)  =  (1.035)  10~7  lb  •  sec/ft2 

u(T  )  =  (6.601)  10~7  lb  •  sec/ft2 
w 

follow.  The  exponent  w  is  then  w  =  0.8611. 


For  this  value  of  and  =  6,  figure  13.8  of  Schlichting  (1968)  gives, 

approximately. 


But 


p6Ud 


~  -  1.16  [w(T6)/(p6U6x)]"1/2 


so  the  parameter  for  laminar  (adiabatic)  flow  is  expressed  in  terms  of  and 
the  Reynolds  number  based  on  x. 


For  example,  if  1  *  1.0  inch  then  the  values  of  corresponding  to  PQ  *  (700, 
1400,  2100)  psia  TQ  -  1100  °R  are  Mt  =  (0.1547,  0.1301,  0.1176).  If  laminar- 
turbulent  transition  is  assumed  to  take  place  within  an  infinitesimally  narrow  x- 
interval  corresponding  to  the  start  of  the  rough  part  of  the  plate,  then  the 
starting  conditions  for  the  solution  of  the  first  order  ordinary  differential 
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equation  for  M^(x)  given  at  the  end  of  section  VI  above  may  be  defined  such  that 
the  starting  at  the  begining  of  the  turbulent  boundary  layer  correspond  to  the 
same  value  of  the  momentum  thickness  as  at  the  end  of  the  laminar  region. 

C.  Specific  Calculations 

The  constants  of  intergration  and  C2  (which  were  introduced  in  section  III  B 
above)  arose  during  the  evaluation  of  a  y  integral.  Accordingly,  they  are  inde¬ 
pendent  of  y.  They  may  however  depend  on  the  parameters  and  3^.  Bradshaw 
(1977)  suggests  the  functions 

ct  =  1.0 

C?  =  5.0  +  95.0  M  2  +  30.7  3  +  226.0  3  2 

1  t  q  q 

A  few  other  constants  not  yet  specified  were  used  in  the  calculations.  They  are 

(i)  The  Karmen  constant:  <  *  0.41 

(ii)  The  ratio  of  specific  heats  for  air:  Y  =  1.4 

(iii)  The  gas  constant  for  air:  R  =  1716.48  (ft.  lbs)/(slug  °R) 

(iv)  The  Coles  wake  parameter  for  flat  plate  flow:  II  *  0.62 

The  derivative  with  respect  to  of  the  function  9(M^)  was  approximated  numeri¬ 
cally  by  a  formula  of  the  form 

dQ  6((1  +  e)MTJ  -  61(1  -  e)MtJ 
___  «  2"T~M 

T  1 

with  e  =  0.05. 

Figures  1  and  2  illustrate  the  effect  of  surface  roughness  height  on  the  distri¬ 
butions  of  skin  friction  and  momentum  thickness,  respectively  for  adiabatic  wall 
conditions.  Figures  3,  4,  and  5  Illustrate  the  effect  of  wall  temperature  on  the 
distributions  of  skin  friction  parameter,  all  for  a  fixed  roughness  height  k  * 
0.04  Inches. 
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The  calculation  results  of  figures  1-5  were  all  based  on  the  same  assumptions 
regarding  tunnel  stagnation  temperature  and  pressure.  These  were  Tq  -  1100  °R  and 
PQ  =  1400  psia,  respectively. 

The  results  agree  with  the  trends  that  one  might  reasonably  expect.  The  only 
thing  that  might  be  surprising  to  readers  not  intimately  familiar  with  high  speed 
boundary  layers  is  the  small  numerical  value  of  the  momentum  thickness.  This  can 
be  rationalized  to  some  extent  by  noting  that  the  absolute  temperature  of  the 
fluid  on  the  wall  is  about  nine  times  that  of  the  fluid  at  the  outer  edge.  Thus 
the  mass  density  of  the  fluid  on  the  wall  is  about  one  ninth  that  of  the  fluid  in 
the  free  scream.  From  Che  definition  of  the  momentum  thickness 


one  sees  that  a  small  value  of  p(y)  in  the  region  where  U/U^  is  intermediate 
between  its  extremes  tends  to  make  the  value  of  8  smaller  than  it  would  be  for 
uniform  density. 

V.  CONCLUSIONS  AND  RECOMMENDATIONS 

Despite  the  crudeness  of  some  of  the  assumptions  employed  in  the  derivation  of  the 
above  algorithm  (such  as,  for  example,  the  neglect  of  energy  losses  due  to  acous¬ 
tic  radiation  and  the  wave  drag  of  individual  roughness  elements)  the  results 
appear  quite  plausible,  and  suggest  that  simple  algorithms  may  be  quite  adequate 
for  the  prediction  of  statistically  two  dimensional  flat  plate  boundary  layers. 

We  feel  that  our  objective  of  furnishing  a  self-contained  derivation  of  an  effi¬ 
cient  and  physically  reasonable  boundary  layer  prediction  method  has  been 
achieved.  Further  work  should,  however,  address  the  neglected  physical  effects 
mentioned  in  the  preceeding  paragraph. 


Norman,  Oklahoma 
June  30,  1984 
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f  1221 < count  > -count*  < 5 ♦ -count > * . 5-1 . 
f  1 1 2  ( cot mt )  =  ( count-2 ♦  >*< count-!  .  >*.5+1 . 
c2 ( mtaul * be tad > -5 ♦ +95 , *mtaul*mtaul 
+betad*  ( 30 . 7+226 .  *betad  > 
mu(  femp>=muref*(  ( temp/tref  >**1 . 5) 

* ( t r  ef +s 1 ) / ( temp+s 1 > 

asi nh<  x>-log<  x+sqrr  ( 1  . +x*x>  > 

ii.fi  me  ( lip  1  us  >  =ii.plus*li  plus/  ( 1  rough*!  rough 

dur  uf  <  !i.p  1  u s )  = .  5*as  1  nn  ( ,  5*l<  f  u nc  (Kpl u s  >  >  //a r  ma > > 

+  ( log(  1  rough ) /li.arman-3 .  >*ii.func(l;.plus> 

/  ( 1 .  +ii.f  unc  ( ii.p  1  us ) ) 
do  3  1  =  1  ?  5  *  2 
z<.  i  +  1  >  — — z  ( i  > 
w  (  i  + 1 ) ~w ( i  ) 


print** "this  program  predicts  the" 
print** "downstream  development  of  a 
print** 'flat  plate  turbulent  boundary' 


-33- 


print*?  7 layer  including  effects  of' 
print*?  compressibil :  ty?  roughness?  ‘ 
prints? 'and  heat  transfer.  to  use  it? ' 
print* •  trie  user  must  specify  whet 
pr  int*?  '  type  of  thermal  boundary  ' 
print*? 'condition  applies7 
tOovtd-l ,  + ,  5*  ( gamma- 1 ,  >  *mde  1 t a*mde 1 t a 
t  •  d  e 1 1 a  — 1 0 / 1 Oo v td 

print*?  type  0  if  the  wall  is  adiabatic?' 

print*? 'type  1  if  the  wall  : n  not' 

pr i nt*  ? ■ ad i aoat ic....' 

read*?  sdwall 

i f ( adwa 1 1 . eq . 0 )  then 

oetsq-0 . 

tw-tde  1  la*  ■  I  ?  + « 5*  ( gamma"  1 «  > *pntur b*mdol  t a 
*moel ta ) /cl 
else 

mt*?  'since  you  have  indicated  that' 
print*? 7 the  wall  is  nonadiahsiio*  a 
print*? 'wall  temperature  is  r ecu i red. 7 
print*? 'type  it  in  (assuming  degrees' 
pr  i  nt*  ?  r  o> »k  i  ne  > « ?  ♦ 
read*?  tw 
end  if 

print*? 'the  current  value  in  storage7 
pr  3  nt* ? 7  for  the  tunne 1  r esevo i r  7 
pr  mt*  ?  7  temperature  is'?t0?7  degrees  rant  me 
print*? 'the  test  section  maeh  number  is' 
pr  ml*  ?mdelta 

prim,*? 'type  in  the  tunnel  resevcir  pressure 
pm  it*  ? 7  (  i  n  ps  i a  > .  *  * 7 
read*  ?  pO 
pO~pO*l 4* ■ 

poo  1  ta~pO*tOovtd**  <  gamma/  ■  i  ■>  -  gamma  >  > 
r bode  i  -pde  1  ta*t-Oovtd/  ( rgas*tO  > 

print*? 'type  in  the  initial  value  of  the' 
pr  i  nt*  ?  7  per  erneter  mta:  '-rode  1  t.;*sqr  t  < cf/2  > 7 
read*  ?  mtau ( 1 ? 

print*? 7 type  in  the  integration  step  sire 
print* ?  m  i nches « ♦ ? 7 
read*  ?  r  i 
i  i~h/12« 

print*? 'type  in  the  number  of  steps' 
print*? 7 over  which  output  data  is  7 
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pr i nt* » ' r  equ i red « ♦ .  ' 
re  sd*  f  n 

print*' 'what  is  the  equivalent  sand' 
print*?  'gram  roughness  he i ght? ( i n 
pr  i  nt*  *  '  i  nodes  >  ? ' 
r  ead* » kr  ough 

ougir"k  r  ou.gh/ 12 . 
steps  “2 

print*' 'what  is  the  length  of' 
print*' 'the  initial  laminar-  run 
pr  1 1  it* »  '  (in  i  nehes )  ? ' 
read*  f  ala": 

;v:iaor-;.-;lam/12, 

uoe  1  ta-mde  1  ta*  sqr  t  ( gamm  a*r  gas*  t-del  t  s 
print*'  'u  oelta-"  »u  delta'  '  tdel  ta='  '■'•del 
print*' '  'xi-am*  '  rhodel-  ' r- node  1 

mom  1  am--  <  58*sqr  t  ( mu  ( toe  1  ta  >  •*  a:  1  am  /  ( ude  1 ta 
*ihodel>) 

print*' '  the  momentum  thickness  a*' 
print*' 'the  end  of  the  laminar  run' 
pr i nt* » ' (call ed  mom 1 am >  i a ' » mom 1 am ' ' f ee 
8  continue 

;>eg i n  v outer)  runge-uvutta  loop 
do  2  1=1 » steps 
count”! < 
mtau2:-mtau  ( i 
theta <  i  >  =0 < 
chmtau-0 , 

•aegiu  calculation  of  the  slope  function  f<mt.a 
which  is  the  right  hand  aide  a*'  the  first 
order  ordinary  differential  equation  to  be 
solved 

I  obetaq  ( i > -0 * 

dtheta-0 . 
do  5  m-~l » 1 ' 2 

m  v$u  1  =mta<  i2*  ( 1 » +esma  1 1  *r e a  1  ( m )  > 

upidel-mdelta*sqrt  ( tdelta/tw  >  ,'mtaul 

1 v i sc= sqr  t ( rgss*t  w/ gamma  >  *mu ( tw > / < pde 1 1 

kp  1  us-t-.r  ough/ 1  vi  sc 

gr  oup!  =  (  gamma."  1 4  >¥mtaul*mtau!: 

group2-sqrt( 4  5*pntur  b*gr  oup 1 > 

i f  < sows 1 1 ♦ ne . 0 )  then 

bo-nag  - ( 'delta* < l .  +  * 5* < gamma- 1 .  > *pnturb 
*moelta*mdel ta)  /tv-cl )  /<pnturb*upi del  '• 


obet-sq  ( i  >  -:obet3q  (  i  )  + .  5*'bet3q  j 

end if  1 

gr oup3=:sqr t  ( c  1  + ,  5*pnt ur b*bet-3q*bet  qq  /dr  oupi  >  I 
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C: 

c  oegin  gqusseqn  qu3dr3fu.ro  to  cqleuiqte 
c  momentum  'thickness 

1  f  < >30 1  p  Is  ♦  I  e « 80 « >  then  ! 

print.)?>  "dubious  deipls.  del  pis-  '  fdeipi:: 

pr  1  nt* ?'m="?m?"  i-'ri 

print*? ' count- " ? count? '  kpius” ' » kpius 

pr  1 1  it*  ?  '  group2- '  ?  group2 

end  if 

ntgrl=0«  1 

••is>moo®c2  <  mtqu  1  ?  be  t  qq  >  + 1  og  <  qw«*  1 1  >  /kqrmqn 
d-dsmoo-dur uf  ( kp  1  us  > 
c  1  nv=  ( 1 ,  - ,  cj/  ( kqr  m-3n*3W3 1 1  >  >  /dsmoo 

_i--0 

eogeb 1 ® .  f  3 1 se * 
yplbot~0 » 
ypl top=25. 
mqxstp=15 
10  continue 

OO  A  I;  -  1  .»  -  | 

y  P 1  - ,  5*  <  z  <  k  >  *  ( yp  1  top-yp  1  bot  >  ; 

+ypltop+yplbot>  ! 

incupi-d*t3nh(ypl*cinv>+(3sinh( «5*ypl/qwqII > 

+co  1  es*  ( 1  *  -cos  <  p  i  *yp  1  /de  1  p  1  s ) ) )  /kqr  rusvi 
up 1 us® sqr t  <  c 1 >  *s i n < gr oup2* i ncupl > /gr  oup2 
+betqq* <  1 .  -cos (group2*i ncupl  >  >/grou.pl 
denom=pntur  b*i  ip  1  us*  ( betqq+  ♦  5*gr  gup  1  *1  ?p  1  us )  +c  1 
ntgrqn®uplus*  ( 1 .  -uplus/upidel  >/denom 
A  ntgr  1  =:ntgr  1  + « 5*  ( y  pi  top-yp  1  bot  >  *u  ( k ) 

*ntgrqn 
yp 1 bot -yp 1 top 
yp] top-yp 1 top*  2 - 
j=j+1 

i  f  (j*  gt .  msp/stp )  ther; 
pr  mt*i "  thetq  integral  not" 
print* ? " complete  qf  ter " » mqxstp • 

"  steps,  quit  calculation" 
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i “steps 
go  to  11 

e 1 sei f < edged l.eq. . true . )  then 
•go  to  11 
end  if 

if  (ypltop.ge.delpls)  thev 
edged 1 t  r  ue « 
y p 1 top-de 1 p 1 s 
end  if 
•So  to  !.  0* 

the  to  1  - 1  v  i  sc*nt  gr  1  f  toe  1  to/  ( tu'Mi  tp  1  da  1  > 
i f  < count . eq ♦ 1 . > then 

the  to  <  i  >  -thetd  <!>-*-•»  5ft het-gl 
end  if 

dtheto-dthet d+  -  5*resl  < m > *thetdl / < esflidi  l*fl'.to<.u 
if  (i.eq. steps)  go  to  2 
f  ~Pito:-»2*mtou2/  ( Aide  l  to*flidel  t s*dthetd  > 

crt.-h*  f  * .  5 

chin  tssu “choit  su +f  1 22 1  <  count )  *er k/3  • 

i f  < co u nt . eq . 4 . >  go  to  ? 

mt du2"mtsu  (  i  >  +f  1 1 2 <  count )  *cr ■  •; 

count -count* 1 . 

go  to  1 

continue 

Kitdu (i+1 ) “iitdu <  i  >  +chrotsu 
continue 

pr  i  nt '  ( 5;-i  *  o  *  1  Op;  j » 3  ,  ?;•/  *  3  *  ?:•••:  '  *  '  flitou '  ? 

' thetd ' f  ' bet sq 

pri  nt  '  ( 1 >  4f  12.9) ' *  <h*resl  <  i  >  *mt3u(  i  >  * 
thets  (  i )  *  obe tsq  (  i  >  *  i  - 1  .♦  steps ) 
if  ( steps . 1 e . 2 >  then 


■f  ■«  r  r  f 


thet-s 


print**'  type  1  if  t 
print** 'is  close  to  mon low.  type  0' 
printer  "if  the  ogreenient  is  poor' 
read** iirotr 

i  f  (  i  trs;tr  .  eq .  0 )  then 
print** 'type  in  3  better  value  of' 
pr i nt* * '  flitsu (I ).,... ' 
read**nitau<  1 ) 

•go  to  ? 
else 

steps  --V 
go  to  8 


J 
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e  iid  if 
end  if 

i  f  ( adua 1 1  ^  no ,  0  :  t;  ioi  i 

wr  i  to  ( 1 » ''  ( l x ,  a .  f  7 ,  1 > 7  > 7  gi  ven  tw-  7  ,  t- 
end  if 

1  r,e  1 1  •  7  ( 1  x  1 a ,  i2 1  a » f  5 .  0 » a »  f  7 ,5)')'  adwa  1 : 
•qbwqi  i .» 7  pO  <  in  psiq)  = 7  •  pO/i  44.  »  7  nitau  <1  n 
r>it  3u  ( 1  > 

1 te ( 1  > "  ( Ixr a»f'->« 4> 7  >  "krough”7  r froudh-fi; 
wr  i  to  (if7  ( Six  r  a  f  1  Ox  f  a  1 7x »  a » 7x » q ) 7  > 7  x 7  f  7  :ut 
f Mot a 7 1  " betqq 7 

write';  i  >  7  ( lx? 4f  12.9 ) 7 )  <h*reql  (  i  > » nitau  ( i ; 

theta ( i > » obetqq ( i > * i =1  * n  > 
printer 'are  you  finished? (type  0' 
print**  7  if  you  are,  type  1  if  you 7 
printer 'want  to  do  another  example 7 
r  ead* » done 

if  < do no. no, 0)  do  to  c 
endf ile(l) 
r  ew  1  nd  ( 1 ) 
end 
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